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ABSTRACT.  We  present  a systematic  approach  to  the  computation  of  exact  nonreflecting 
boundary  conditions  for  the  wave  equation.  In  both  two  and  three  dimensions,  the  crit- 
ical step  in  our  analysis  involves  convolution  with  the  inverse  Laplace  transform  of  the 
logarithmic  derivative  of  a Hankel  function.  The  main  technical  result  in  this  paper  is 
that  the  logarithmic  derivative  of  the  Hankel  function  H{vx\z)  of  real  order  v can  be  ap- 
proximated in  the  upper  half  z -plane  with  relative  error  £ by  a rational  function  of  degree 
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term  coming  from  the  calculation  of  a spherical  harmonic  transform  at  each  time  step.  In 
short,  nonreflecting  boundary  conditions  can  be  imposed  to  any  desired  accuracy,  at  a cost 
dominated  by  the  interior  grid  work,  which  scales  like  A2  in  two  dimensions  and  A3  in 
three  dimensions. 
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NONREFLECTING  BOUNDARY  KERNELS 


1.  Introduction 

A longstanding  practical  issue  in  numerical  wave  propagation  and  scattering  problems 
concerns  the  reduction  of  an  unbounded  domain  to  a bounded  domain  by  the  imposition 
of  nonreflecting  boundary  conditions  at  an  artificial  boundary.  We  restrict  our  attention  to 
“time-domain”  calculations,  for  which  it  is  well-known  that  the  exact  nonreflecting  condi- 
tions are  global  in  both  space  and  time.  While  the  problem  has  been  widely  studied  (see 
Givoli  [ 1 ] for  an  overview),  the  boundary  conditions  used  in  practice  typically  introduce 
serious  numerical  artifacts.  The  two  most  common  approaches  are  based  on  the  construc- 
tion of  local  differential  boundary  conditions  [2,  3]  or  absorbing  regions  [4,  5],  but  neither 
provides  a clear  sequence  of  approximations  which  converge  to  the  exact,  nonlocal  con- 
ditions. Recently,  Sofronov  [6]  and,  independently,  Grote  and  Keller  [7]  have  developed 
and  implemented  an  integrodifferential  approach  for  three-dimensional  calculations  using  a 
spherical  boundary  and  have  demonstrated  that  high  accuracy  can  be  achieved  at  reasonable 
cost.  In  their  schemes,  the  work  is  of  the  same  order  as  the  explicit  finite  difference  or  finite 
element  calculation  in  the  interior  of  the  domain.  For  TV2  points  on  the  spherical  boundary, 
0(TVJ)  work  is  required.  Hagstrom  and  Hariharan  [8]  have  shown  that  these  conditions 
can  be  effectively  implemented  using  only  local  operators,  but  at  the  cost  of  introducing  a 
large  number  of  auxiliary  functions  at  the  boundary.  A somewhat  more  general,  but  closely 
related,  integral  formulation  is  introduced  in  [9,  10].  The  fundamental  analytical  tool  in  the 
latter  papers  is  what  we  refer  to  as  the  “nonreflecting  boundary  kernel”  which  is  the  inverse 
Laplace  transform  of  the  logarithmic  derivative  of  a Hankel  function. 

In  this  paper,  we  prove  that  the  logarithmic  derivative  of  a Hankel  function  can  be  ap- 
proximated as  a ratio  of  polynomials  of  modest  degree,  so  that  its  inverse  Laplace  transform 
can  be  expressed  as  a sum  of  exponentials.  Our  analytical  approach  combines  an  exten- 
sion of  the  Mittag-Leffler  theorem  with  the  approximation  techniques  of  the  fast  multipole 
method.  In  particular,  Theorem  4. 1 presents  an  exact  representation  of  the  logarithmic  de- 
rivative as  a sum  of  poles  plus  a continuous  density  on  the  branch  cut.  Theorem  4.6,  which 
is  preceded  by  several  technical  lemmas,  presents  a reduced,  approximate  representation. 

Using  this  approach,  the  cost  of  computing  the  nonreflecting  boundary  condition  is  com- 
parable to  that  of  a fast  Fourier  or  spherical  harmonic  transform.  For  two-dimensional 
problems,  O (TV  log  TV  log  i)  work  is  required  at  each  time  step,  where  TV  is  the  number  of 
points  used  in  the  discretization  of  a cylindrical  (circular)  boundary.  In  three  dimensions, 
the  cost  is  proportional  to  TV 2 log2  TV  + TV2  log  TV  log  j,  for  a spherical  boundary  with  TV2 
points.  The  first  term  comes  from  the  calculation  of  the  spherical  harmonic  transform  using 
the  fast  algorithm  of  [1 1,  12]. 

Other  authors,  including  Nedelec  [13]  and  Cruz  and  Sesma  [14],  have  studied  the  log- 
arithmic derivative  of  the  Hankel  function,  based  on  a variety  of  techniques.  In  this  paper 
we  present  a sum-of-poles  representation  for  the  logarithmic  derivative  of  a Hankel  func- 
tion of  real  order  v bounded  away  from  zero  with  accuracy  s for  argument,  z,  satisfying 
Im(z)  > 0.  The  number  of  poles  is  bounded  by  0[  log  |v|  -log  j -blog2  |v|-Hv|-1  log2  j).  A 
similar  representation  for  v = 0 is  also  derived  which  is  valid  for  Im(z)  >rj>  0 requiring 
0(  log  | • log  j + log  j • log  log  \ + log  i • log  log  i)  poles. 
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In  Section  2,  we  introduce  nonreflecting  boundary  kernels.  In  Section  3 we  collect  back- 
ground material  in  a form  convenient  for  the  subsequent  development.  Section  4 contains 
the  analytical  and  approximate  treatment  of  the  logarithmic  derivative,  while  a procedure 
for  computing  these  representations  is  presented  in  Section  5.  The  results  of  our  numerical 
computations  are  contained  in  Section  6,  and  we  present  our  conclusions  in  Section  7. 


2.  Nonreflecting  Boundary  Kernels 


Let  us  first  consider  the  wave  equation 

utt=c2\'2u  (1) 

in  a two-dimensional  annular  domain  po  < p < p\.  The  general  solution  can  be  expressed 
as 

oo 

u(p,<p,t)=  ein<t>  C~x[an{s)Kn(ps/c)  +bn(s)  In(ps /c)](t),  (2) 

n=— oo 


where  Kn  and  In  are  modified  Bessel  functions  (see,  for  example,  [15]  9.6), 

Kn(z)  = ^in+lH^(zeni/2),  In(z ) = rnJn(zem/1),  -tt  < argz  < y,  (3) 

the  coefficients  an  and  bn  are  arbitrary  functions  analytic  in  the  right  half  plane,  C denotes 
the  Laplace  transform 


-j£ 


oo 


C[f](s)=  I e~st  f (t)  dt, 


(4) 


and  C 1 denotes  the  inverse  Laplace  transform 


1 fl0C 

£-1[g](0  = z—  / estg(s)ds. 
2 m J — i oo 


(5) 


Likewise,  for  the  wave  equation  in  a three-dimensional  domain  ro  < r < r\ , the  general 
solution  can  be  expressed  as 


oo  n 


u(r,  4>,  6,  t)  = ^ Ynm((p,d)jC  1 


n=— oo  m=—n 


■ Kn+,_(rs/c) 

&nm  (,S  ) / -j— 

y/rs/c 
In+i{rs/c)-\ 


(6) 


4"  bnm  (s) 


y/rs/c 


if). 


If  we  imagine  that  p = p\  (or  r = r\)  is  to  be  used  as  a nonreflecting  boundary,  then  we 
can  assume  there  are  no  sources  in  the  exterior  region  and  the  coefficients  bn  is)  (or  bnm  (s)) 
are  zero.  Let  us  now  denote  by  u„(p,  t ) the  function  satisfying 


C[un](p,  s ) = a„(s)  Kn(ps/c). 


(7) 
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Then 


^[^««](/>.  s)  = an(s ) • S-  ■ K'n{ps/c ) 

= C[un]tP,s).(S-K>s/c)\ 
\c  Kn(ps/c)) 


(8) 


so  that 


—un(p,  t)  = un(p , t)*£ 
dp 


s_  K'n(ps/c ) 
c Kn(ps/c)_ 


(0, 


where  * denotes  Laplace  convolution 


/' 


(f*g)(0=  / - z)dz. 


(9) 


(10) 


The  convolution  kernel  in  (9)  is  a generalized  function.  Its  singular  part  is  easily  removed, 
however,  by  subtracting  the  first  two  terms  of  the  asymptotic  expansion 


5 K'„(ps/c)  _ s 1 , , 

c Kn{ps/c)  c 2p  v >' 


00. 


(11) 


From  the  assumption  un(p,  t)  = 0 for  t < 0 and  standard  properties  of  the  Laplace  trans- 
form we  obtain  the  boundary  condition 


a 13  1 

— u„(p , t)  + un(p , t ) + —un(p,  t ) = 

3p  c 3r  2p 


fo, 

Jo 


(t  - t)  un(p,  z)dz , 


where 


a„(0  = £ 


-l 


£ J_  £^r'(p5/c) 


(0, 


(12) 


(13) 


_c  2 p c Kn(ps Jc) _ 

which  we  impose  at  p = p\ . 

Remark.  The  solution  to  the  wave  equation  in  physical  space  is  recovered  on  the  nonre- 
flecting boundary  from  un  by  Fourier  transformation: 

N/2-\ 

u(p\,<j),t)=  ^ Un(p\,t)ein<p,  (14) 

n=-N/2 

assuming  N points  are  used  in  the  discretization. 

The  analogous  boundary  condition  in  three  dimensions  is  expressed  in  terms  of  the  func- 
tions unm(r,  t)  satisfying 

Kn+l  irs/c) 

C[unm]{r,s)  = anm(s) -== — . (15) 

y/rs/c 

After  some  algebraic  manipulation,  assuming  unm  (r,  t)  = 0 for  t < 0,  we  have 

a 

Jr 


i a i r 1 

unm{r,t)  + -—unm(r,t)  + -unm{r,t)  = / a)„(t  - z)  unm(r,  z)  dz, 
cot  r J o 


(16) 
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where 


ton  (t)  — C' 


-1 


5 1 5^+.("A)-| 

-4 1 

c 2 r c Kn+\_(rs/c) _ 


(0, 


(17) 


which  we  impose  at  r = r\ . 

Note  that  the  boundary  conditions  (12)  and  (16)  are  exact  but  nonlocal,  since  they  rely  on 
a Fourier  (or  spherical  harmonic)  transformation  in  space,  and  are  history  dependent.  The 
form  of  the  history  is  simple,  however,  and  expressed,  for  each  separate  mode,  in  terms  of 
a convolution  kernel  which  is  the  inverse  Laplace  transform  of  a function  defined  in  terms 
of  the  logarithmic  derivative  of  a modified  Bessel  function 


-7-I0  gKv(z) 

dz 


K(z) 

Kv{z)' 


(18) 


Remark.  In  three  dimensions,  the  required  logarithmic  derivative  of  Kn+\{z)  is  a ratio 
of  polynomials,  so  that  one  can  recast  the  boundary  condition  in  terms  of  a differential 
operator  of  order  n.  The  resulting  expression  would  be  equivalent  to  those  derived  by 
Sofronov  [6]  and  Grote  and  Keller  [7]. 


The  remainder  of  this  paper  is  devoted  to  the  approximation  of  the  logarithmic  deriva- 
tives (1 8)  as  a ratio  of  polynomials  of  degree  O (log  v),  from  which  the  convolution  kernels 
on  and  <jon  can  be  expressed  as  a sum  of  decaying  exponentials.  This  representation  allows 
for  the  recursive  evaluation  of  the  integral  operators  in  (12)  and  (16),  using  only  O(\o>gn  ) 
work  per  time  step  (see  [16]).  We  note  that,  by  Parseval’s  equality,  the  Z2  error  resulting 
from  convolution  with  an  approximate  kernel  is  sharply  bounded  by  the  Zoo  error  in  the 
approximation  to  the  kernel’s  transform.  Precisely,  approximating  the  kernel  B(t ) by  the 
kernel  A(t)  we  find 


A * u — B * u |L  = II  Au  — Bu  |L  < sup 


A - B 


seiR  1 5 1 

A-B I 


= sup  — — 

s M B 


Bu 


B *u  |L 


(19) 


where  we  assume  that  A,  B,  and  u are  all  regular  for  Re (5)  > 0.  For  finite  times  we  may 
let  5 have  a positive  real  part,  r] : 


\\A*u- B*u\LA(!T)<e,<T  sup  l|£*»lLz(0.n-  (20) 

serj+iR  | D | 

We  therefore  concentrate  our  theoretical  developments  on  Zoo  approximations.  For  ease  of 
computation,  however,  we  compute  our  rational  representations  by  least  squares  methods. 
These  do  generally  lead  to  small  relative  errors  in  the  maximum  norm,  as  will  be  shown. 
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Since  Hankel  functions  are  more  commonly  used  in  the  special  function  literature,  we 
will  write  the  logarithmic  derivatives  as 


d log  Kv(z)  = ~~  log  {z  e7tl/2)  = i Hv 


(1  )'/_7ri 


(zeni/2) 


dz 


Hi])(ze^/2) 


(21) 


We  are,  then,  interested  in  approximating  logarithmic  derivative  of  the  Hankel  function  on 
and  above  the  real  axis. 


3.  Mathematical  Preliminaries 


In  this  section  we  collect  several  well-known  facts  concerning  the  Bessel  equation,  the 
logarithmic  derivative  of  the  Hankel  function,  and  pole  expansions,  in  a form  that  will  be 
useful  in  the  subsequent  analytical  development. 


3.1.  Bessel’s  Equation.  Bessel’s  differential  equation 

,2 


d2u  1 du  / y^\ 


dz2 


(22) 


for  y g R,  has  linearly  independent  solutions  and  //u(2),  known  as  Hankel’s  functions. 
These  can  be  expressed  by  the  formulae 

/*«(*)  = (23) 


r(1)  ^ J.v{z)  - e-™*  Jv{z) 


i sin(v7r) 

where  the  Bessel  function  of  the  first  kind  is  defined  by 

- 00  (-z2/ 4)k 


i sin(v7r) 


J"{Z)  ^ §*!r(v + * + !)■ 


(24) 


The  expressions  in  (23)  are  replaced  by  their  limiting  values  for  integer  values  of  v.  (See, 
for  example,  [15,  9.1])  For  general  v,  the  functions  //v(1)  and  //y2)  have  a branch  point  at 
z = 0 and  it  is  customary  to  place  the  corresponding  branch  cut  on  the  negative  real  axis 
and  impose  the  restriction  — tt  < argz  < i r.  We  shall  find  it  more  convenient,  however,  to 
place  the  branch  cut  on  the  negative  imaginary  axis,  with  the  restriction 

71  3n 

--  < argz  < — . (25) 

Hankel’s  functions  have  especially  simple  asymptotic  properties.  In  particular  (see,  for 
example,  [17,  7.4.1]), 


■> 7ZZ - 


Hl'hz)  ~ ( 


7 Tz' 


Ak(v) 

yk  ’ 

(26) 

k= 0 

k=0 

(27) 
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as  z — > oo,  with  — it  + 8 < argz  <2 jt  — 8,  where 

(4v2  - 12)(4v2  - 32)  • • • (4v2  - (2£  - l)2) 


^A-(v)  = 


A'!  8* 


and  the  branch  of  the  square  root  is  determined  by 
Finally  we  note  the  symmetry 


zl/2  _ e{\og\z\+i  argz)/2 


(28) 

(29) 


= (30) 

We  also  make  use  of  the  modified  Bessel  functions,  Kv(z ) and  Iv(z).  These  are  linearly 
independent  solutions  of  the  equation  obtained  from  (22)  by  the  transformation  z ->  iz. 
Their  Wronskian  satisfies 

Kv(z)I'v(z)  — K'v(z)Iv(z)  = z~l.  (31) 

Moreover  we  have  for  positive  r [ 1 8] 


H^\re~in/1)  = —e~vni/2(ev7tiKv{r)  + 7ziIv(r)). 
ni 


(32) 


Asymptotic  expansions  of  Kv(r ) and  Iv(r)  for  r small  and  large  are  also  known  [15,  9.6, 
9.7].  For  real  r and  v > 0 we  have 


Kv(r) 

IV(T) 

Kv{r) 

Ur) 


Y-  log-,  v = 0. 


in©’- 


V(v+  1)V2- 
'-jz 


1 


2jzr 


er , 


0, 

0, 

oo, 

00. 


(33) 

(34) 

(35) 

(36) 


Here  y = 0.5772  ...  is  the  Euler  constant. 

Finally,  we  note  the  uniform  expansions  of  Bessel  functions  for  v -*  oo  given  in  [15] 
For  Hankel  function  and  derivative  we  have 

4£  \ i/4Ai(e27r//3v2/3^) 


H^(yz)  ~ 2^— (y~Z~T2  ) 


,1/3 


H<iy(vz) 


4e~2ni^  / 4£  \-l/4Ai/(e2jn73v2//3£) 


.2/3 


as  v —>  oo,  where  we  restrict  z to  |arg(z)|  < jt/2  and  define 


2 3/2  . 1 + Vl  - Z2  r 

-c  7“  = log v 1 — z2. 

3 z 


(37) 

(38) 

(39) 
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Here,  Ai(t)  denotes  the  Airy  function  [15,  10.4].  Note  that  £ = 0 when  z — 1.  Large  v 
approximations  of  the  modified  Bessel  functions  for  real  arguments,  r,  are  given  by 


Kv{yr)  ~ J 


where 


zv  (!  +r2 )V4’ 

</>(r)  = log 


/v(vr) 


s/2nv  (i  +/-2)1/4’ 


oo. 


1 + Vl  + r2 


+ vT+r2. 


(40) 


(41) 


3.2.  Hankel  Function  Logarithmic  Derivative.  We  denote  the  logarithmic  derivative  of 

H{VX)  by  G„, 


Gv(z ) = log  H^\z)  = - (~} 

dz  Hf\z) 


(42) 


The  following  lemma  states  a few  fundamental  facts  about  Gv  that  we  will  use  below. 
Lemma  3.1.  The  function  Gv(z),for  v e R satisfies  the  formulae 

G-V(z ) = Gv(z),  (43) 


Gv(zeni)  = Gv(z)en\ 


Tt  Tt 

< argz  < — , 


(44) 


Gv(z) 


z->0,  (45) 


2 ° ~ 2 

where  z —\z  \ e~l  arg~  A r/ze  complex  conjugate  of  z.  Asymptotic  approximations  to  Gv  are 

[\og(ze~ni/1 /2)  + y)  1 z~]  + 0(z),  v = 0, 

— |y | z— 1 + 0(z2^~l),  0 < \v\  < 1, 

-|u|z_1  + O(zlogz),  |v|  = 1, 

— |v|  z-1  + 0(z),  \v\  > 1, 

where  y is  the  Euler  constant, 

k= 0 * k= 0 

where  A/fi >)  is  defined  in  (28),  and 

2e~ni^  / 4£  \-i/2Ai'(e27n/3v2/3£) 


oo, 


(46) 


Gv(vz) 


V1/3z 


&)■ 


Ai(e27n/3v2/3£)  ’ 


oo, 


where  £ is  defined  in  (39).  Furthermore,  the  function  uv  defined  by 

uv(z)  = z Gv(z) 

satisfies  the  recurrence 


(47) 


(48) 


uv(z)  = 


V - 1 - U — i (z) 


— V. 


(49) 
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Im(z) 


FIGURE  1 . Curve  z(£ ) defined  by  (39)  near  which  the  scaled  zeros  of  //„(1  * 
lie  (see  Lemma  3.2).  The  branch  cut  of  is  chosen  (25)  on  the  negative 
imaginary  axis. 


Proof.  Equations  (43)  and  (44)  and  asymptotic  expansion  (45)  follow  immediately  from 
the  definitions  (23)  through  (25)  of  Jv  and  Hy  . The  asymptotic  expansion  (46)  follows 
from  (26)  and  (27),  while  (47)  is  a consequence  of  (37)  and  (38).  The  recurrence  (49)  is 
derived  from  standard  Bessel  recurrences  [15,  9.1.27],  □ 

The  zeros  of  H 1 ) (z)  are  well  characterized  [15,  18];  they  lie  in  the  lower  half  z-plane  near 
the  curve  shown  in  Figure  1 obtained  by  transformation  [19]  of  Bessel’s  equation.  In  terms 
of  the  asymptotic  approximation  (37),  this  curve  corresponds  to  negative,  real  arguments 
of  the  Airy  function. 

Lemma  3.2.  The  zeros  hVi\ , hvg,  • ■ ■ of  Hy\z)  in  the  sector  —tz/2  < argz  < 0 are  given 
by  the  asymptotic  expansion 


K«~vz(in)  + 0(v-'),  „ = 1,.'..,LM/2  + l/4j, 

uniformly  in  n,  where  is  defined  by  the  equation 

Kn=e-^i\-^  an. 


(50) 


(51) 


z(£)  is  obtained  from  inverting  (39),  and  an  is  the  nth  negative  zero  of  Airy  function  Ai. 
The  zeros  in  the  sector  tz  < argz  < 3tz/2  are  given  by  —hv  \ , —hvg, In  particular, 


~ i;  + e-2,n/3(v/2)1/3(-a,). 


(52) 


where  —a\  = 2.338 

3.3.  Pole  Expansions.  A set  of  poles  in  a finite  region  defines  a function  that  is  smooth 
away  from  the  region,  with  the  smoothness  increasing  as  the  distance  increases.  This  fact 
leads  to  the  following  approximation  related  to  the  fast  multipole  method  [20,  21]. 
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Lemma  3.3.  Suppose  that  q 1 , . . . , qn  are  complex  numbers  and  z\ , . . . 
numbers  with  \zj  | < 1 for  j — 1 , . . . , n.  The  function 


n 


/(z)  = T, 

7 = 1 


can  be  approximated  for  Re(z)  = a > 1 by  the  m pole  expansion 


772  -1 


gm(z)  = 

7=0 


Yj 

z — aP  ' 


where  co  = e2jzi/!m  is  a root  of  unity  and  Yj  Is  defined  by 


J 772  — 1 72 

Yj  = — Y'  0)~Jl  y ^ qkZk,  j = o, . . . , m - 1. 

m f-r'  “ 


/=0  /t=l 

77?c  error  o/ the  approximation  is  bounded  by 


2{cr  + 1) 


|/w-*,«|<(a-_1)(a_1)2if«i. 


where 


F(z)  = £fiL, 

jZ 1 Z Z7 


Proof  We  use  the  geometric  series  summation 


^ = £ 


W—  1 y772  | 


+ 


z — o 


2t=0 


^Tr+l  Z — V 


to  obtain 


772—  1 i 72  772  — 1 

/(z)  - g„  (z)  = £ gj  z/'  - Y,  vj  ) 

k=0  j= 1 7=0 


, zn  are  complex 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 


(59) 
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All  m terms  of  the  first  summation  vanish,  due  to  the  combination  of  (55)  and  the  equality 
°^k  — m <$A-0-  For  the  error  term  we  obtain 


and 


n ^izim 


j=\ z 


n I <tizjm 


EV7 -7  <y ^ i7_i7_  < _ 

7 7 _•  < ^ I 7 7 _•  1 7 1 | •,  / 

r = l |1  ZjjZ 


7 = 1 


‘7 


I " 

Izl 


1^7 


7: 


<l^±L±(l^Ll}M<^±LlRe(Y 

\z\  (a  - l)2  1 + a 2 


\qj\ 


{a  - l)2  |z|  V " 1 -Zj/z 


) 


< 


I?/ 1 


a2  + 1 


{a  - \)2  \ - Zj 


EE 


a2  + 1 
(a-  l)2 


|F(z)|, 


1 y Y] 

\m~l  V 

T Yj 

1 Z — (x)j 

1 4-i,  z — co-l 

7=0 


7=0 


= |gm(E)|. 


(60) 


(61) 


Moreover,  repeating  the  computations  of  (60),  we  find 

\m\  < a'  + 1 |F(z) |.  (62) 

(«  - lr 

Now  the  combination  of  (59)  through  (62)  and  the  triangle  inequality  gives  (56).  □ 


Inequality  (56)  remains  valid  if  we  assume  instead  that  |zy  | < b and  Re(z)  = ab  > b,  for 
arbitrary  b £ R,  b > 0;  this  fact  leads  to  the  next  two  results  whose  proofs  mimic  that  of 
Lemma  3.3  and  are  omitted. 


Lemma  3.4.  Suppose  n,  p are  positive  integers,  q \, ...  , qn  are  complex  numbers,  and 
z\ , . . . , zn  are  complex  numbers  contained  in  disks  D\ , . . . , Dp  of  radii  r\ , . . . , rp,  cen- 
tered at  c\, . . . ,cp,  respectively.  The  function 


n 


f(z)  = J2 

7 = 1 


(63) 


can  be  approximated  for  z satisfying  Re(z  — a)  > art  > rl  for  i = 1 , . . . , p by  the  m ■ p 
pole  expansion 


p w-i  v. . 

gm(z)  = TEZ 77y 

frtjro  z~(ci+  n or) 


(64) 


where  gij  Is  defined  by 


Yij 


■iS*""  E.  "-(A2)' 


/=0  zkeDi\Ui-\ 


i = 1 , . . 

2=0,.. 


• ,p, 

• , m - 1, 


(65) 
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with  Ui  = U j<jDj.  The  error  of  the  approximation  is  bounded  by 


where 


\f{z)  - gm(z)\  < 


2(a2  + 1)  \F(z)\ 
(am  - l)(fl  - l)2’ 


(66) 


n 


n 2)  = £ 

j= l 


(67) 


Lemma  3.5.  Suppose  that  the  discrete  poles  of  Lemma  3.4  are  replaced  with  a density  q 
defined  on  a curve  C with  C C Up  = D\  U • • • U Dp,  specifically 


/CO 


-I 


q(0 


dt;. 


<cz-S 

which  is  finite  for  z outside  Up,  and  that  gm  is  defined  by  (64)  with  yij  defined  by 

m—\  n . 

i = \, ...  ,p, 

j =0, ...  ,m  - l. 


(68) 


Yu  = <»~il  f M—1^ 


(69) 


with  Ui  = U j<iDj.  Then  the  bound  (66)  holds  as  before. 


Lemma  3.3  enables  us  to  approximate,  with  exponential  convergence,  a function  defined 
as  a sum  of  poles.  The  fundamental  assumption  is  that  the  region  of  interest  be  “separated” 
from  the  pole  locations.  The  notion  of  separation  is  effectively  relaxed  by  covering  the  pole 
locations  with  disks  of  varying  size  in  an  adaptive  manner.  In  Lemmas  3.4  and  3.5,  we  use 
this  approach  to  derive  our  principal  analytical  result. 


4.  Rational  Approximation  of  the  Logarithmic  Derivative 


The  logarithmic  derivative  of  the  Hankel  function  Gv(z)  approaches  a constant  as  z 
oo  and  is  regular  for  finite  z e C,  except  at  z = 0,  which  is  a branch  point,  and  at  the  zeros 
of  //y(1)(z),  all  simple.  We  can  therefore  develop  a representation  for  Gv  analogous  to  that 
of  the  Mittag-Leffier  theorem:  the  only  addition  is  due  to  the  branch  cut  on  the  negative 
imaginary  axis.  It  will  be  convenient  to  work  with  uv(z),  for  which  approximations  to  be 
introduced  have  simple  error  bounds. 


Theorem  4.1.  The  function  uv(z)  = z Gv(z),  where  Gv  is  defined  for  v e R by  (42)  with 
the  branch  cut  defined  by  (25),  is  given  by  the  formula 


Nv 


uv(z)  = iz 


i hv.n  l 

7 ^ 7 ~ hV'„  I*’ 


n= 1 


i r°° 

ni  Jo 


Im (uv{re  ni^ 2)) 


ir  + z 


dr 


(70) 


for  z e C not  in  {0,  hv  \.  hvp, . . . , hv  ^v}  and  not  on  the  negative  imaginary  axis.  Here 
hv,i,  hv,  2, . . . ,hVNv  denote  the  zeros  of  H^\z),  which  number  Nv. 
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Figure  2.  Integration  contour  Cm,  with  inner  circle  radius  1 /m  and  outer 
radius  m + 1 . 


Proof.  The  case  of  the  spherical  Hankel  function,  where  v = k + 1/2  for  k e Z,  is  simple 
and  we  consider  it  first.  Here  uv(z)  is  a ratio  of  polynomials  in  iz  with  real  coefficients, 
which  is  clear  from  the  observation  that  u 1/2(2)  = iz  — 1/2  in  combination  with  the  recur- 
rence (49).  Hence 


Ny 


Uv(z)  = p(z ) + av'"  , 

^ 7 - hVt„ 


n= 1 


where  p is  a polynomial  and  ofv,„  is  the  residue  of  uv  at  hVM, 

(Xv.n  = lim  (z  hvn ) uv(z)  — hvn 


(71) 


(72) 


by  THopitaTs  rule.  We  see  from  (46)  that 

uv(z)  ~ iz  — - -(-  0(z_1),  z — >■  00,  (73) 

whence 

p(z)  = iz-j.  (74) 

Noting  that  uv(iy)  e E for  y € 1R,  and  combining  (71),  (72),  and  (74),  we  obtain  (70). 

We  now  consider  the  case  v ^ k + 1/2,  k e Z,  for  which  the  origin  is  a branch  point. 
For  m = 1, 2, . . . , we  define  Cm  to  be  the  simple  closed  curve,  shown  in  Figure  2,  which 
proceeds  counterclockwise  along  the  circle  ] of  radius  m + 1 centered  at  the  origin  from 
argz  = — 7r/2  to  3jt/2,  to  the  vertical  segment  z = rei7t1^2,  r e [1/m,  m + 1],  to  the 
circle  C„7  of  radius  1/m  centered  at  the  origin  from  argz  = 3jt/2  to  —n/2,  to  the  vertical 
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Figure  3.  Plot  of  Re(wv(re  ni^2)),  containing  the  zero  crossing,  and 
Im(wy(re_jri//2)),  for  v = 2 and  r e [0,  3]. 


segment  z = re  ni^2,  back  to  the  first  circle.  Since  none  of  the  zeros  of  H^V)  lies  on  the 
imaginary  axis,  Cm  encloses  them  all  if  m is  sufficiently  large.  For  such  m,  and  z e C 
inside  Cm  with  Hv(z)  ^ 0,  the  residue  theorem  gives 


i r uv(£) 

2 7Ti  Jcm  s -Z 


Nv 

dt;  = uv(z)  + 

n= 1 


h 


hv,n 

v,n  2 


(75) 


We  now  consider  the  separate  pieces  of  the  contour  Cm.  For  the  circles  C^1}  and  C^  \ we 
use  the  asymptotic  expansion  (73)  about  infinity  and  (45)  about  the  origin  to  obtain 


1 

lim  

w-+oo  2rci 


dt;  = iz  - -, 


lim  

2ni 


d;  = o. 


(76) 


Now  exploiting  the  symmetry  uv{re2ni^2)  — uv(re~ni /2)  from  (44)  for  the  vertical  seg- 
ments, we  obtain 


1 

lim  — 

m->oc  2 Tli 


f Uv(t 

' Jcm  C - 


uv (K)  i 

-dt;  =iz--  + 
z I In 


Lf 

m Jo 


oo 


2i\m{uv(re  711 Z2)) 


( re  ni /2  — z) 

which,  when  combined  with  (75),  yields  (70)  and  the  theorem. 


-ni/1  dr , 


(77) 

□ 


The  primary  aim  of  this  paper  is  to  reduce  the  summation  and  integral  of  (70)  to  a similar 
summation  involving  dramatically  fewer  terms.  To  do  so,  we  restrict  z to  the  upper  half- 
plane and  settle  for  an  approximation.  Such  a representation  is  possible,  for  the  poles  of  uv 
(zeros  of  //v(1))  lie  entirely  in  the  lower  half-plane  and  do  not  cluster  near  the  real  axis.  We 
first  examine  the  behavior  of  uv  on  the  negative  imaginary  axis. 

The  qualitative  behavior  of  uv  on  the  branch  cut  is  illustrated  by  the  case  of  v = 2,  shown 
in  Figure  3.  The  plot  changes  little  with  changing  v,  except  for  the  sign  of  Im(wv(z))  and 
the  sharpness  of  its  extremum. 
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Lemma  4.2.  For  v e M,  v k + 1/2,  k e 1,  the  function  uv(re~ni/1)  is  infinitely 
differentiable  on  r e (0,  oo)  and  has  imaginary  part  satisfying  the  following  formulae: 


\m(uv{re  711 /2)) 


IT  COS (VTZ) 

cos2(vir)K2(r)  + (irlv(r)  + sin(y^)ATv(r))2 


# 0, 


(78) 


JTCOS(VTT)  2|U| 


Im  (uv(re~m/2)) 

4M-1  r(|y|)2 

Im (uv(re~7:i^2))  ~2cos (v7r)re~2r, 


IT 


(log(r/2)  + y)2  +7T2' 


v = 0, 
v 0, 


Im (uv(re  ni,2j) 
where  (p  is  defined  in  (41). 


COS  (V7T)fr2  + V2 


cosh  (2v  0(r/|y|))  + sin(|  y |tt)  ’ 


r —>  0, 

r — ► oo, 
M ->  00, 


(79) 

(80) 
(81) 


Proof  Infinite  differentiability  of  uv(z)  follows  from  the  observation  that  H{hX) (z)  0 on 

the  negative  imaginary  axis.  To  derive  (78)  we  recall  (32)  to  obtain 


Im (uv(re-*i/2))  = r7r  cosjv  n ){Kv(r)rv(r)  - A»/„(r)) 

cos2(vir)K2(r)  + {jt  I v{r)  + sin(y7r)AT(r))“ 


(82) 


then  apply  (31).  The  remaining  formulas  follow  from  the  asymptotic  forms  of  Kv(r)  and 
Iv(r ) for  small  and  large  r,  and  the  uniform  large  v expansions  given  in  (33)  through  (36) 
and  (40).  Here  we  use  the  symmetry  u-v  = uv.  Note  that  (79)  is  valid  for  r/ |v|  — > 0.  The 
approximation  (81)  is  nonuniform  for  v % 2k  — 1 /2  and  irlv(r)  + sm(vir)Kv(r)  % 0.  □ 


Lemma  4.3.  Given  vo  > 0 there  exist  constants  Co  and  c\  such  that,  for  all  v e K,  |y|  > vq, 
v k -f  1 /2,  k e Z,  and  all  z satisfying  Im  (z)  > 0,  the  function 


f(z)  = 


Im (uv(re  Z2)) 

dr 

ir  + z 


(83) 


satisfies  the  bounds 


cp 

1 + M/M 


< 1/00 1 < 


Cl 

1 + \z\/\v\- 


(84) 


Moreover,  there  exists  8 > 0 such  that  for  all  v e E,  \v\  > yo,  and  e with  0 < e < 1/2, 
f (z)  admits  an  approximation  g(z)  that  is  a sum  of  d < <5-(l  + |y|  1 log(l/s))  -log(l/f) 
poles,  with 


|/(z)-g(z)|  <£•  |/(z)|, 


(85) 


provided  Im(z)  > 0. 
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Proof.  We  assume  v ^ k + 1/2  for  integral  k and  begin  by  changing  variables,  r = |v| w, 
so  that 


lm(uv(\v\we  m/2)) 
iw  + z/\v\ 


du)  = 


ixz(w)dw. 


(86) 


From  the  nonvanishing  of  fiz  and  its  asymptotic  behavior  in  w,  it  is  clear  that  (84)  holds 
for  |u|  € (no,  ui)  and  any  fixed  v\  > v0-  Using  (81)  for  |u|  large  but  bounded  away  from 
2 k — 1/2  for  integral  k,  an  application  of  Watson’s  lemma  to  (83)  focuses  on  the  unique 
positive  zero,  w*,  of  cp  defined  in  (41).  As  the  derivative  of  this  function  is  positive,  we 
conclude 


a cos(u7r) 
iw*  + z/\v\ ’ 


(87) 


where  a is  a function  of  w*,  so  that  (84)  clearly  holds.  However,  as  v 2k  — 1/2,  the 
denominator  on  the  right-hand  side  of  (81)  may  nearly  vanish  at  w*  and  the  expansion 
loses  its  uniformity.  Setting  cos(v7r)  = r]  in  these  cases,  we  see  that  the  denominator  has  a 
minimum  which  is  bounded  below  by  0(r]2).  Hence  in  an  0(|v|-1)  neighborhood  of  the 
minimum  which  includes  w*  we  have 


J A h(w) 


rj\v \y/\  + ( w *)2  1 

I ds 

iW*+z/\v\  J-y/\v\  VT  + f2V2S2 


(88) 


which  by  the  change  of  variables  5 = rjz/ |v|  is  seen  to  satisfy  the  upper  bound  in  (84) 
uniformly  in  r].  As  the  rest  of  the  integral  is  small,  the  upper  bound  holds. 

We  now  move  on  to  the  approximation.  For  a positive  integer  m and  a positive  number 
wo,  we  define  intervals  7o  = (0,  wo),  Ij  = ( 2J~lwo , 2-^wo)  for  j = 1, . . . , m,  and  Im+\  = 
(2mwo,  oo).  Now 


f(z)  = fo(z)  + f\  (z)  + fl(z). 


(89) 


where  fo,  f\,  and  fz  are  defined  by  the  formulae 


fo(z)  = 


/xz(w)dw. 


/iz(w)dw. 


(90) 


We  will  now  choose  wo  and  m so  that  fo  and  fz  can  be  ignored  and  then  use  Lemma  3.5 
to  approximate  f\ . Using  (79)  and  (81)  and  taking  wq  sufficiently  small  we  have,  for  some 
constant  cz  independent  of  v, 


l/o(2)  I < 


_C2\v\_(3e  2lv\  fW0 
1 + \z\/\v\  4 Jo 


w2M~]dw  < 


c 2 


1 + \z\ 


•(?“*) 


2\v\ 


(91) 


Hence,  a choice  of 


WW),  € ->  0, 


WQ  = 0(8 


(92) 
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suffices  to  guarantee 


l/oOOl  < |l/(*)l,  (93) 

in  the  closed  upper  half-plane.  Now  using  (80)  and  (81)  and  assuming  m sufficiently  large 
we  have,  for  some  constant  c3  independent  of  v. 


l/2(2)|  < 

From  (92),  choosing 


c3|y| 


1 + \z\/\v\  J2mwo 


roo 

J2m  Wt 


we 


- i>  u; 


dw  < 


C2 


1 + \z\/\v 


■2mw0e-^imw\  (94) 


m > mo  + m\—  log-, 

M £ 


(95) 


for  appropriate  m o and  m i independent  of  v and  e leads  to 

I/2C0I  < |l/(z)|.  (96) 

Finally,  we  apply  Lemma  3.5  to  the  approximation  of  f\  The  error  involves  the  function 
F\  — f |Im(w„)|/0>  + z)dr,  but  we  note  that  |Fi  | = \f\ |.  Using  p poles  for  each  j we 
produce  a p ■ m-pole  approximation  g(z)  with  an  error  estimate,  again  for  Im(z)  > 0,  given 
by 


1/iU)  ~g(z)\  < 


37TTl/i«l- 


(97) 


A choice  of 

p = 0(  logf),  (98) 

enforces 

|/i(z)  -g(z)l  < ||/(z)|.  (99) 

By  combining  (93),  (96),  (99),  and  the  triangle  inequality,  we  obtain  (85)  with  the  number 
of  poles,  d = p ■ m,  satisfying  the  stated  bound.  □ 


The  case  v = 0 requires  special  treatment.  First,  the  direct  application  of  the  preceding 
arguments  leads  to  a significantly  larger  upper  bound  on  the  number  of  poles.  Secondly, 
we  note  that  «o(0)  = 0,  so  that  relative  error  bounds  near  z = 0 require  a vanishing 
absolute  error.  Finally,  the  lack  of  regularity  of  uq(z)  at  z = 0 precludes  uniform  rational 
approximation,  as  discussed  in  [9].  Therefore,  we  relax  the  condition  Im(z)  > 0 to  Im(z)  > 
rj  > 0.  By  (20)  this  will  lead  to  good  approximate  convolutions  for  times  T < . 


Lemma  4.4.  There  exists  8 > 0 such  that  for  all  s,  0 < s < 1 /2  and  rj,  0 < rj  < 1 /2, 
the  function  f{z)  = uq(z)  — iz  + 1/2  admits  an  approximation  g(z)  that  is  a sum  of 
d < 8 • ( log( I/77)  + log log(l / sf)  ■ log(l/f)  poles,  with 

\f(z)-g(z)\<£-\f(z)\. 


provided  Im(z)  > 77. 


(100) 
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Proof.  Note  that  since  w0(O  has  no  poles,  / (z)  is  given  by  (83),  and  satisfies  (84).  Define 
intervals 

/ = ((27_1  ~ 0*7.  (2y  ~ 0*7)  for  j = 1 m,  /m+1  = (( 2m  - 1)77,  00). 

Now 


/(z)  = /i(z)  + /2(z), 

where  /1 , and  /2  are  defined  by  the  formulae 


/i(0 


_ f Im(up(rg  ^/2)) 


y=i  ^ 


ir  + z 


dr,  f2(z ) 


= / 

J Im+\ 


lm.{uQ{re  nil2)) 


ir  + z 


001) 


dr.  (102) 


We  will  now  choose  m so  that  /2  can  be  ignored  and  then  use  Lemma  3.5  to  approximate 
/1 . Using  (80)  and  assuming  m sufficiently  large  we  have,  for  some  constant  c, 

»oo 


1/2(01  < 


1 + \z\ 


L 


— 2r  7 ^ 

rg  aid  < 

1 + \z\ 


2m-lne-2^ 


(103) 


Hence,  choosing 


m > m0(log(l/O  + loglog(l/g)),  (104) 

for  appropriate  mo  independent  of  rj  and  s leads  to 

1/2(01  < |l/(OI-  (105) 

Finally,  we  apply  Lemma  3.5  to  the  approximation  of  f\ . Using  p poles  for  each  j we 
produce  a p ■ m- pole  approximation  g(z)  with  an  error  estimate  for  Im(z)  > r]  given  by 

|/l (z)-g(z) | < ^ryl/^OI-  (106) 


A choice  of 

p = 0(  log-),  (107) 

£ 

enforces 

l/i  (O  -g(OI  < | l/(OI-  (108) 

By  (105),  (108),  and  the  triangle  inequality,  (100)  is  achieved  with  the  number  of  poles, 
d = p ■ m,  satisfying  the  stated  bound.  □ 

We  now  consider  the  contribution  of  the  poles. 

Lemma  4.5.  There  exist  constants  Co,  C\,  8 > 0 such  that  for  all  v,  e e M with  2 < |v| 
and  0 < e < 1/2  the  function 


h(z)  = 

n=  1 


(109) 
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where  hvj, . . . ,hv^v  are  the  roots  of  H^  \ satisfies  the  inequalities 

C\\v\  C2I 

,,,  ■ < M*)\  < T— 7 777  7,  (1  10) 

1 -Hz|/|y|  \ + \z\/\v\ 

and  admits  an  approximation  g(z ) that  is  a sum  of d < <5  • log  |v»|  - log(l/£)  poles,  with 

\h(z)  — g(z)\  < s • \h(z)\  (111) 

provided  Im(z)  > 0. 


Proof  The  curve  C defined  in  Lemma  3.2,  near  which  hVt i/|v|, . . . ,hv^v/\v\  lie,  is  con- 
tained in  disks  separated  from  the  real  axis.  If  we  denote  the  disk  of  radius  r centered  at  c 
by  D{r,  c ),  then  the  disks 

{£>(— Im(z),  z)|  z € C,  | argz  — jt/2|  = jt/2  + n/2n,  « = 1,2, ...},  (112) 

for  example,  contain  C\{+1,  —1}.  From  (52),  the  root  hvj  closest  to  the  real  axis  satisfies 

arg^v,i  ~ ' (113) 

hence  is  contained  in  a disk  of  (112)  with  n ~ log2  (24/33'  a\)  1 1 v |2/3),  and  all 

of  the  roots  are  contained  in  O (log  |v|)  of  the  disks.  Now  applying  Lemma  3.4  we  obtain 
(111)  with  \h  \ replaced  by  \H\  — \ J2  fiv.n  |/(z  — hVM)\.  To  obtain  the  upper  bound  in  ( 1 1 0) 
for  both  h and  H we  note  first  that  it  is  trivial  except  for  |z/v|  ~ 1.  A detailed  analysis  of 
the  roots  as  described  by  Lemma  3.2  shows  that 

|lm(A„j)|  > c//3|u|'/3.  (114) 

Hence,  for  |z/v|  % 1, 


E 


hvj 


z - hv , 


< C|v|2/3^y_2/3  < 3C|y|. 
j= 1 


(115) 


The  lower  bound  in  (1 10)  is  again  obvious  except  for  |z/v|  ~ 1.  Then,  however,  we  note 
that 

h(z)  = uv(z)  — iz  + 1/2  — /(z).  (116) 


Since,  from  (47),  \uv{z)\  = 0(|y|2/3)  for  |z/v|  ~ 1 and  |/(z)|  = 0(1)  by  (84)  the  right- 
hand  side  is  dominated  by  —iz  and  \h(z)\  = (9(|y|).  □ 


The  combination  of  Theorem  4. 1 and  Lemmas  4.3  and  4.5  suffices  to  prove  our  principal 
analytical  result. 

Theorem  4.6.  Given  vo  > 0 there  exists  8 > 0 such  that  for  all  v e R,  |v|  > vo,  and 
0 < s < 1/2  there  exists  d with 

d < <5( log | y | • log(l/£)  4-log2  |y|  + |y|_1  log2(l/^)). 


(117) 
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and  complex  numbers  a\, . . . ,otd  and  ft, . . . , ft/,  depending  on  v and  e,  such  that  the 
function 


Uv,e(z)  = iz-^  + J2~- 

z n= 1 z 

approximates  uv(z)  with  the  bound 

| Uv(z)  - t/v>e(z)|  < £ • |«v(z)|, 
provided  that  Im(z)  > 0.  Furthermore 

/ r00’  X 1/2 

/ r°°  X 1/2 

|mv(*)  - £fte(.x)|  dxj  <£-lj  \uv(x)  - ix  + \/2^  dx\ 

Proof.  We  first  note  the  lower  bound 


(118) 


(119) 


(120) 


\uv(z)  -iz+  1/2|  > C|V|  . (121) 

l + |z|/|u| 

For  y > 0 the  function  is  nonvanishing  and  has  the  correct  asymptotic  behavior,  so  we  need 
only  consider  the  case  of  |v|  large.  The  result  then  follows  from  (47).  This  proves  (120) 
and  (119)  with  uv  replaced  by  uv  — iz  + 1/2  on  the  right-hand  side.  From  (47)  we  have 

| uv(z)  - iz  + 1 /2 1 < cjy|1/3|u„(z)|,  (122) 

so  that  the  final  result  follows  from  the  scaling  s — ^ |y|  1/3f.  □ 


The  number  of  poles  in  (1 17)  required  to  approximate  uv(z)  to  a tolerance  £,  depends  on 
both  £ and  v.  The  asymptotic  dependence  on  £ is  proportional  to  |y|-1  log2(l/f).  We  will 
see  in  the  numerical  examples,  however,  that  this  term  is  important  only  for  small  |v|;  oth- 
erwise the  dominant  term  is  the  first,  for  an  asymptotic  dependence  of  0[  log  \ v\  -log(l/£)). 
As  we  generally  have  £ M-1  in  practice,  the  term  log2  |v|  is  of  less  importance. 

Similarly,  Lemma  4.4  leads  to  the  following  theorem  for  v = 0. 

Theorem  4.7.  There  exists  8 > 0 such  that  for  all  £,  0 < £ < 1/2  and  rj,  0 < r]  < 1 /2 
there  exists  d < <5-(  log(l/?7)-log(l/£)-)-loglog(l/£)-t-loglog(l  /tj))  and  complex  numbers 
a\ , . . . , oid  and  ft , . . . , ft,  depending  on  rj  and  £,  such  that  the  function 

1 d 

y0,,(2)  = /z--  + T— V U23) 

2 Z Pn 


approximates  uo(z)  with  the  bound 

| Uq(z)  - Uo,e(z)\  < £ • \uq(z)\. 


(124) 
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provided  that  Im(z)  > rj.  Furthermore 


oo 


2 \1/2 

|w0U  + iri)  - Uo,e(x  + irj)\~dx  ) 


oo 


< £ 


oo 


— oo 


\ 1/2 

|tfvO  + irj)  - ix  + Tj  + \/l\2dx  ) . (125) 


Proof.  Again  we  already  have  (124)  with  uq{z)  - iz  + 1/2  on  the  right-hand  side.  By  (45) 
we  find 


| uo(z)  - iz  + 1/2|  < c\og(\/r])\u0(z)\.  (126) 

The  theorem  follows  from  the  scaling  e — > log-1  (1  /t])s.  □ 

As  we  must  take  rj  = T-1,  we  see  that  the  number  of  poles  required  may  grow  like 
log(l  Je)  • log  T + log  T ■ log  log  T.  However,  this  is  only  for  the  mode  n = 0 in  the  two 
dimensional  case.  In  short,  the  T dependence  is  insignificant  in  practice. 


5.  Computation  of  the  Rational  Representations 


Analytical  error  bound  estimates  developed  in  the  previous  sections  are  based  on  max- 
imum norm  errors  as  in  (19)  and  (20).  In  numerical  computation  it  is  often  convenient, 
however,  to  obtain  least  squares  solutions.  Our  method  of  computing  a rational  function 
Uv<£  that  satisfies  ( 1 1 9)  is  to  enforce  (120).  An  alternative  approach  would  be  to  use  rational 
Chebyshev  approximation  as  developed  by  Trefethen  and  Gutknecht  [22,  23,  24]. 

In  the  numerical  computations,  we  work  with 

uv(z)  = uv(z)  — iz  + 1/2  (127) 


and  its  sum-of-poles  approximation  UVt£(z)  = UVi£(z)  — iz  + 1/2.  In  particular,  we  have 
the  nonlinear  least  squares  problem 


min 

p.Q 


Pix) 

Q(x) 


UV(x) 


2 

dx. 


(128) 


for  P,  Q polynomials  with  deg(P)  -I-  1 = deg(0  = d.  Problem  (128)  is  not  only  nonlin- 
ear, but  also  very  poorly  conditioned  when  P,  Q are  represented  in  terms  of  their  monomial 
coefficients.  We  apply  two  tactics  for  coping  with  these  difficulties:  linearization  and  or- 
thogonalization. 

We  linearize  the  problem  by  starting  with  a good  estimate  of  Q , and  updating  P,  Q 
iteratively.  In  particular,  we  solve  the  linear  least  squares  problem 


min 

p(i+\)Q{i+\) 


p«+'Hx) 

WF) 


e(i+1)w  - „ , 

UV(X) 

Qil)(x) 


2 

dx, 


(129) 


where  the  integral  is  replaced  by  a quadrature.  The  initial  values  P(0),  Q(0)  are  obtained  by 
exploiting  the  asymptotic  expansion  (46)  and  the  recurrence  (49).  We  find  that  two  to  three 
iterations  of  (129)  generally  suffice. 
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The  quadrature  for  (129)  is  derived  by  first  changing  variables, 

/oo  pTi/ 2 m 

f{x)dx  = / /(tan#)  sec2  6 d9  % iy,-  /(tan#,)  sec2#/,  (130) 

-oo  J—n/2  z_i 

where  #i , . . . ,6m  and  iui  . . . ,wm  denote  appropriate  quadrature  nodes  and  weights.  The 
transformed  integrand  is  periodic  on  the  interval  \—n/2,  n/2],  so  the  trapezoidal  rule  (or 
midpoint  rule)  is  an  obvious  candidate.  The  integrand  is  infinitely  continously  differen- 
tiable, except  at  # = 0,  where  its  regularity  is  of  order  2|v|.  For  |v|  >8  (say),  the 
trapezoidal  rule  delivers  at  least  16th-order  convergence  and  is  very  effective.  For  small 
|v|,  however,  a quadrature  that  adjusts  for  the  complicated  singularity  at  # = 0 is  needed. 
Here  we  can  successively  subdivide  the  interval  near  the  singularity,  applying  high-order 
quadratures  on  each  subinterval  (see,  for  example,  [25]).  We,  instead,  apply  quadratures 
designed  to  handle  a variety  of  singularities  efficiently  [26]. 

The  quadrature  discretization  of  (129)  cannot  be  solved  as  a least-squares  problem  by 
standard  techniques,  due  to  its  extremely  poor  conditioning.  We  avoid  forming  the  corre- 
sponding matrix;  rather  we  solve  the  least  squares  problem  by  Gram-Schmidt  orthogonal- 
ization.  The  2d  + 1 functions 


U y,  1,  XUV,  X , 

. . . , xd  ]uv,  xd  x duv 

(131) 

are  orthogonalized  under  the  real  inner  product 

</.£>(=  / 

r°°  Re(/(x)g(x))  ^ 

-oo  \Qii)(x)\2 

(132) 

to  obtain  the  orthogonal  functions 

Uv(x), 

n = 1, 

gn(x)  = • 

1, 

n = 2, 

(133) 

, x v^rnin{4,« 

Xgn- 2{X)  - 2^=1 

~]}  cnj  gn-j(x),  n = 3, ...  ,2d+  1, 

where 

{xgn— 2i  gn—j)i 

CnJ  ~ i \ 

(gn—j ? gn—j)i 

n = 3, . . . , 2d  + 1, 

7 = 1,---  , min{4,  n - 1). 

(134) 

Now 

g2d+\  = - 

-P°'+1)  + uv  Qii+]\ 

(135) 

so  and  (2(z+1)  are  computed  from  the  recurrence  coefficients  cnj  by  splitting  (133) 

into  even  and  odd-numbered  parts. 

For  some  applications,  including  nonreflecting  boundary  kernels,  it  is  convenient  to  rep- 
resent P / Q as  a sum  of  poles, 


P(z) 

Q(z) 


E 


<*n 

Z ~ Pn‘ 


(136) 


NONREFLECTING  BOUNDARY  KERNELS 


23 


We  compute  p\,  ...  , fid  (zeros  of  Q)  by  Newton  iteration  with  zero  suppression  (see,  for 
example,  [27])  by  the  formula 


Aa+1)  = Au) 


(137) 


where  , . . . , are  the  previously  computed  zeros  of  Q.  Then  cx\, . . . , ad  are  com- 
puted by  the  formula  an  = P(Pn)/Q'(Pn)-  The  derivative  Q'(z ) is  obtained  by  differenti- 
ating the  recurrence  (133). 


6.  Numerical  Results 


We  have  implemented  the  algorithm  described  in  Section  5 to  compute  the  representa- 
tions of  on  and  con  through  their  Laplace  transforms.  Recall  that  for  the  cylinder  kernels, 
on,  we  have  v — n while  for  the  sphere  kernels,  con,  we  have  v = n + 1/2.  Table  1 presents 
the  sizes  of  the  representations  for  s = 10-6,  10-8,  and  10“ 15  in  (120).  For  the  cylinder 
kernels,  which  are  affected  by  the  branch  cut,  the  number  of  poles  for  small  n is  higher 
than  for  the  sphere  kernels.  This  discrepancy,  however,  rapidly  vanishes  as  n increases  and 
the  asymptotic  performance  ensues.  The  log(l/£)  dependence  of  the  number  of  poles  for 
n > 1 0 is  clear. 

For  s = 10-8  we  have  also  computed  the  maximum  norm  relative  errors  which  appear 
in  (19)  by  sampling  on  a fine  mesh.  For  the  cylinder  kernel  with  n = 0,  we  expect  an  <9(1 ) 
error  in  a small  interval  about  the  origin  due  to  (79).  However,  errors  of  less  than  e are 
achieved  for  |s|  > 5 x 10~7.  This  implies  a similar  accuracy  in  the  approximation  of  the 
convolution  for  times  of  order  106.  For  all  other  cases  the  maximum  norm  relative  errors 
are  of  order  e. 

Finally,  Table  2 presents  poles  and  coefficients  for  the  cylinder  kernels  for  n = 1 , . . . ,4 
and  e = 10-6,  to  allow  comparison  by  a reader  interested  in  repeating  our  calculations. 
Note  that  the  pole  locations  are  written  in  terms  of  5 = z/i.  Extensive  tables  will  be  made 
available  on  the  Web  at  http : / /math . nist . gov/mcsd/Staf  f /BAlpert. 

Remark.  Our  approximate  representation  of  the  nonreflecting  boundary  kernel  could  be 
used  to  reduce  the  cost  of  the  method  introduced  by  Grote  and  Keller  [7].  The  differential 
operators  of  degree  n obtained  in  their  derivation  need  only  be  replaced  by  the  correspond- 
ing differential  operators  of  degree  log  n for  any  specified  accuracy.  It  is  interesting  to  note 
that  in  the  two-dimensional  case,  where  the  approach  of  [7]  does  not  apply,  the  analysis 
described  above  can  be  used  to  derive  an  integrodifferential  formulation  in  the  same  spirit. 


7.  Summary 


In  this  paper  we  have  introduced  new  representations  for  the  logarithmic  derivative  of  a 
Hankel  function  of  real  order,  that  scale  in  size  as  the  logarithm  of  the  order.  An  algorithm 
to  compute  the  representations  was  presented  and  our  numerical  results  demonstrate  that 
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Table  1 . Number  d of  poles  to  represent  the  Laplace  transform  of  nonre- 
flecting boundary  kernels  on  and  con , for  various  values  of  s. 


£ = 

10~6 

On 

C On 

n 

d 

n 

d 

o 

26 

£ = 

10~15 

1 

9 

On 

C On 

2 

6 

n 

d 

n 

d 

3-6 

5 

0-5 

n 

1 

41 

7-8 

6 

6-8 

6 

2 

24 

9-  12 

7 

9-  12 

7 

3 

18 

13-  19 

8 

13-  19 

8 

4 

15 

20-31 

9 

20-31 

9 

5 

14 

32-51 

10 

32-51 

10 

6 

13 

52-86 

11 

52-86 

11 

7-12 

12 

87-  147 

12 

87-  147 

12 

13-  14 

13 

0-  13 

n 

148-227 

13 

148-228 

13 

15-  16 

14 

14-  15 

14 

228-401 

14 

229-402 

14 

17-18 

15 

16-18 

15 

402-  728 

15 

403-  728 

15 

19—22 

16 

19-21 

16 

729-1024 

16 

729—1024 

16 

23-26 

17 

22-25 

17 

1 8 

ic  'Jf) 

1 8 

£ = 

10~8 

A / J I 

32-37 

1 o 

19 

A\J — jU 

31-36 

1 o 

19 

on 

C On 

38-45 

20 

37-44 

20 

n 

d 

n 

d 

46-  54 

21 

45-53 

21 

0 

44 

55-65 

22 

54-  65 

22 

1 

15 

66-  79 

23 

66-  79 

23 

2 

9 

80-97 

24 

80-96 

24 

3-8 

7 

0-7 

n 

98-118 

25 

97-118 

25 

9-  10 

8 

8-  10 

8 

119-  145 

26 

119-  144 

26 

11-  14 

9 

11-  14 

9 

146-  177 

27 

145-  176 

27 

15-20 

10 

15-  19 

10 

178-216 

28 

177-216 

28 

21-28 

11 

20-28 

11 

217-265 

29 

217-264 

29 

29-41 

12 

29-40 

12 

266-  324 

30 

265-  324 

30 

42-58 

13 

41-57 

13 

325-  397 

31 

325-  396 

31 

59-  84 

14 

58-83 

14 

398-486 

32 

397-485 

32 

85-123 

15 

84-123 

15 

487-  595 

33 

486-  594 

33 

124-  183 

16 

124-  183 

16 

596-  728 

34 

595-  727 

34 

184-275 

17 

184-  275 

17 

729-  890 

35 

728-  890 

35 

276-418 

18 

276—418 

18 

891-1024 

36 

891-1024 

36 

419-638 

19 

419-637 

19 

639-971 

20 

638-971 

20 

972-1024 

21 

972-1024 

21 

the  new  representations  are  modest  in  size  for  orders  and  accuracies  likely  to  be  of  practical 
interest. 
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TABLE  2.  Laplace  transform  of  cylinder  kernel  on  defined  in  (13),  approx- 
imated as  a sum  of  d poles,  for  n = 1, . . . , 4 and  s = 10-6. 


Pole  Coefficient 

Pole  Location 

n 

d 

Re 

Im 

Re 

Im 

1 

9 

-0.426478£  - 02 

0.000000£  + 00 

— 0.368403£  +01 

0.000000£  + 00 

-0.416255£  -01 

0.000000£  + 00 

— 0.205860£  +01 

0.000000£  + 00 

-0.122665£  +00 

0.000000£  + 00 

— 0.1 18994£  +01 

0.000000£  + 00 

-0.143704£  + 00 

0.000000£  + 00 

— 0.717570£  +00 

0.000000£  + 00 

-0.530662 E - 01 

0.000000£  + 00 

— 0.423506£  + 00 

0.000000£  + 00 

—0.863872 E - 02 

0.000000£  + 00 

—0.2231 1 1£  + 00 

0.000000£  + 00 

—0.961472 E - 03 

0.000000£  + 00 

— 0.103710£  +00 

0.000000£  + 00 

— 0.721 548£  -04 

0. 000000 £ + 00 

— 0.409342£  - 01 

0.000000£  + 00 

-0.250102£  - 05 

0.000000£  + 00 

— 0. 1 17 156£  — 01 

0.000000£  + 00 

2 

6 

0.218164£  -01 

0.000000£  + 00 

— 0.333263£  + 01 

- 0.000000£  + 00 

0.860648£  + 00 

0.000000£  + 00 

— 0.162945£  +01 

0.000000£  + 00 

-0.138934£  +01 

0.162069£  +00 

-0.125843£  +01 

0.412637£  + 00 

-0.138934£  +01 

— 0.162069£  +00 

-0.125843£  + 01 

-0.412637£  +00 

0.209905 E - 01 

0.000000£  + 00 

-0.612710£  +00 

0.000000£  + 00 

0.232032 E - 03 

0.000000£  + 00 

— 0.240327£  + 00 

0.000000£  + 00 

3 

5 

-0.179277£  +00 

0.000000£  + 00 

-0.309775£  +01 

0.000000£  + 00 

-0.168335£  +01 

0.1291 1 1£  + 01 

— 0.167998£  + 01 

0.130784£  +01 

-0.168335£  + 01 

—0.1291 1 1£  +01 

-0.167998£  +01 

— 0.130784£  + 01 

-0.816322£  +00 

0.000000£  + 00 

-0.187260£  +01 

0.000000£  + 00 

— 0.126962£  — 01 

0.000000£  + 00 

— 0.950854£  + 00 

0.000000£  + 00 

4 

5 

-0.197725£  +01 

0.220886£  +01 

— 0.197861£  +01 

0.220444£  +01 

-0.197725£  +01 

— 0.220886£  +01 

— 0. 197861  £ + 01 

— 0.220444£  + 01 

-0.219247£  +01 

0.216535£  +01 

— 0.282304£  + 01 

0.382237£  + 00 

-0.219247£  +01 

— 0.216535£  + 01 

-0.282304£  +01 

— 0.382237£  + 00 

0.464435£  + 00 

0.000000£  + 00 

-0.201 159£  +01 

0.000000£  + 00 

The  present  motivation  for  this  work  is  the  numerical  modeling  of  nonreflecting  bound- 
aries for  the  wave  equation,  discussed  briefly  here  and  in  more  detail  in  [16].  Maxwell’s 
equations  are  also  susceptible  to  similar  treatment  as  outlined  in  [28].  The  new  represen- 
tations enable  the  application  of  the  exact  nonreflecting  boundary  conditions,  which  are 
global  in  space  and  time,  to  be  computationally  effective. 
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